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We consider the genesis and dynamics of interfacial instability in gas-liquid flows, 
using as a model the two-dimensional channel flow of a thin falling film sheared by 
counter-current gas. The methodology is linear stability theory (Orr-Sommerfeld 
analysis) together with direct numerical simulation of the two-phase flow in the case 
of nonlinear disturbances. We investigate the influence of three main flow parame¬ 
ters (density contrast between liquid and gas, film thickness, pressure drop applied 
to drive the gas stream) on the interfacial dynamics. Energy budget analyses based 
on the Orr-Sommerfeld theory reveal various coexisting unstable modes (interfacial, 
shear, internal) in the case of high density contrasts, which results in mode coales¬ 
cence and mode competition, but only one dynamically relevant unstable internal 
mode for low density contrast. The same linear stability approach provides a quanti¬ 
tative prediction for the onset of (partial) liquid flow reversal in terms of the gas and 
liquid flow rates. A study of absolute and convective instability for low density con¬ 
trast shows that the system is absolutely unstable for all but two narrow regions of 
the investigated parameter space. Direct numerical simulations of the same system 
(low density contrast) show that linear theory holds up remarkably well upon the 
onset of large-amplitude waves as well as the existence of weakly nonlinear waves. In 
comparison, for high density contrasts corresponding more closely to an air-water- 
type system, although the linear stability theory is successful at determining the 
most-dominant features in the interfacial wave dynamics at early-to-intermediate 
times, the short waves selected by the linear theory undergo secondary instability 
and the wave train is no longer regular but rather exhibits chaotic dynamics and 
eventually, wave overturning. 
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I. INTRODUCTION 

Vertical two-phase flows of thin liquid Aims sheared by a counter-current gas are prototyp¬ 
ical for many technical applications, such as absorption and distillation (using structured 
packing), evaporation and condensation. In these applications, not only mass and heat 
transfer but also operational limits are closely linked to the prevailing hydrodynamics. The 
flow in the two phases, in turn, is largely determined by the interactions between gas and 
liquid at its interface. Although gas-sheared liquid Aims have been part of active research 
on fundamental and practical level for several decades, the rich interfacial dynamics are still 
not fully understood. These rich dynamics emerge as the liquid interface becomes unstable, 
leading to the development of waves. Depending on the flow rates, these waves can form 
highly complex structures, which may give rise to wave break-up, ligament formation and 
droplet entrainment. 

One of the first to investigate a channel flow with two superimposed fluid layers from a 
theoretical point was Yih [1], who used asymptotic expansion to solve the Orr-Sommerfeld 
eigenvalue problem associated with the temporal linear stability in the long wavelength limit 
for equal densities and layer thicknesses. He found that viscosity stratification alone can 
cause interfacial instability at arbitrarily small Reynolds numbers, which is also referred to 
as Yih mechanism. Yiantsios [2] extended the linear stability analysis by accounting for short 
waves as well as effects due to surface tension and gravity. They observed that, depending 
on the choice of parameters, the flow is receptive to a short-wave instability at low Reynolds 
numbers and, moreover, to a shear mode instability (Tollmien-Schlichting mechanism) for 
sufficiently large Reynolds numbers. To classify the various types of instabilities arising 
in parallel two-phase flow, Boomkamp et al. [3] analysed by which mechanism energy is 
transferred from the primary flow to growing disturbances, thereby verifying that both the 
Yih and shear mode are important routes to interfacial instability. In fact, a combination of 
these two mechanisms, also referred to as internal mode, represents a further possible route. 
An established method to obtain the full spectrum of eigenvalues in the Orr-Sommerfeld 
problem is the Chebyshev collocation method |1] . This method can also be used to determine 
whether the interface is convectively or absolutely unstable 13 E], a characteristic that is of 
importance for the operation of processes utilizing shear flows. An extensive review on the 
analysis of spatially developing flows can be found in Reference [7] . 

Linear stability analysis is an effective and proven technique to gain valuable insight in 
the genesis of interfacial instability. However, its validity is limited to disturbances with 
infinitesimal amplitude. As these disturbances grow, nonlinear effects become important 
and have to be taken into account. Due to the complexity of nonlinear stability only few 
general theories exist [8]. Nonetheless, a variety of modelling techniques have been proposed 
to describe the development of the interfacial waves up to finite amplitude. Some of these 
techniques impose a priori assumptions about the wave dynamics under investigation, like 
long-wave or lubrication approximation, resulting in model equations such as the Kuramoto- 
Sivashinsky equation or depth-averaged integral equations. A large number of studies on the 
nonlinear dynamics of interfacial flows are based on these kind of models iHIS]. However, 
given that their range of applicability is generally not known in advance, they may at 
times be prone to ambiguous results PI- This is especially the case for fairly “punishing” 
flow regimes, involving large pressure fluctuations and potentially large-amplitude waves, 
for which there is a major necessity to gain fundamental understanding. On the other 
hand, weakly nonlinear theories based on either the Stuart-Landau or the Ginzburg-Landau 
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equations dispense with assumptions that cannot be conhrmed a priori and are capable 
of matching experimental observations [TH - flG] . Such “pure” weakly nonlinear theories are 
therefore appropriate to complement direct numerical simulations (DNS) of the full Navier- 
Stokes equations, which, in turn, are guided by basic Orr-Sommerfeld (OS) analysis, to 
study interfacial instability in a rigorous manner. 

The aforementioned techniques have helped to shed light on the genesis and development 
of interfacial instability in shear flows and a good understanding of the mechanisms at play 
has been gained. However, a substantial amount of the available literature is dedicated 
to horizontal flows or flows down an inclined plane. Flow dynamics specihc to a vertical 
conhguration have not received as much attention even though the same methods are ap¬ 
plicable. Phenomena related to this conhguration, such as partial liquid how reversal due 
to counter-current gas (hooding), inhuence the design, optimization as well as operation 
of technical processes to a great extent and have hitherto mainly been investigated exper¬ 
imentally. We therefore try to help elucidate the rich dynamics of vertical counter-current 
gas-liquid hows from a more theoretical standpoint using the (semi-)analytical methods and 
simulation techniques described above. 

In this work we consider the laminar-laminar case only, although the ideas and results 
contained herein can be extended to turbulent gas streams. The most accurate methodology 
appropriate for such an enhanced level of complexity is full-scale DNS, which is the target of 
future work. However, short of full-scale DNS, a quasi-laminar model may be assumed for 
the linear stability of the two-phase how mm or, alternatively, a weighted residual integral 
boundary layer method (WRIBL) model [181 - 120] . which also models nonlinear interfacial 
waves. These approaches both have their own advantages and shortcomings, and the aim 
of future work will be to conhrm these reduced-dimensional modelling approaches with 
evidence from accurate DNS, towards which the present work is an initial contribution. 

This work is organized as follows. We give a description of the investigated problem and 
outline the employed methods in § [^ Results regarding temporal stability of the system are 
discussed in l^lHIl Spatio-temporal behaviour with respect to absolute/convective instability 


of the linearized dynamics is presented in § |IV[ Section [V| discusses nonlinear wave dynamics. 
Finally, we give concluding remarks in ^ |VI 


II. PROBLEM DESCRIPTION AND COMPUTATIONAL METHODS 

In this analysis we consider the dynamics of a gas-liquid flow in a vertical channel, de¬ 
scribed schematically in Fig. [Tj The two continuous phases are separated by an, initially, 
flat interface. A pressure gradient dp/dx > 0 in vertical direction counteracts gravity. 
We investigate cases in which the balance between gravity and pressure gradient gives rise 
to counter-current flow, with gas flowing in the direction of decreasing pressure and liq¬ 
uid flowing in the direction of gravity. Figure also shows the development of a linear 
small-amplitude wave at the interface. Typically, the evolution of interfacial waves depends 
sensitively on the details of the mean flow. 

Both fluid layers exhibit steady, spatially uniform, laminar and incompressible flow along 
the vertical rectangular channel. To describe this two-dimensional flow, we use a Cartesian 
coordinate system, {x,z), in which the flat interface is located at z = 0 and the conhning 
channel walls are located at z = —cIl and at z = do, respectively. Within these boundaries, 
the fully developed liquid and gas layer occupy the regions —di < z < 0 and 0 < 2 ; < do, 
respectively. 
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FIG. 1. Schematic representation of the undisturbed base flow. Both fluids are assumed laminar. 
The dashed line indicates the development of a wave at the interface. 


A. Base flowf and linear stability analysis 


With the above mentioned conditions, the Navier-Stokes equations describing the fluid 
flow in both phases reduce to standard balances between pressure, viscous and gravitational 
forces, which are subject to no-slip condition at the channel walls, 2 : = —di and = da, as 
well as continuity of tangential stress at the interface, z = 0. To write the governing equa¬ 
tions in nondimensional form we introduce the following dimensionless variables (without 
tildes) and scalings: 


X = Hx, U = VpU, PgV^ = H — , Tint = PcVlint, t = —t, (1) 

where x = {x, z) and u = {u, w) are the coordinate and velocity vector, H is the channel 
height, Vp is an inertial pressure scale, dp/dx is a positive pressure gradient. Tint is the 
interfacial shear stress, 14,is the gas side interfacial friction velocity and t denotes time. 
Further, the following dimensionless parameters arise: 


Pl Pl ^ dj 

m= —, r = —, dj = —, 

Pg Pg H 


Rcp = 


PGVpH 

Pg 


RCg — 


pcy/^H 


Re-r = 


Pg 


PG^*,intH 

Pg 


We = 


PgVp^H 

7 


( 2 ) 


Here, pj is the dynamic viscosity and pj the density of the respective phase (j = L,G), 
whereas 6j is the relative thickness of the respective fluid layers. The Reynolds numbers 
Rcp, Rcg and Rcr, in turn, relate to the applied pressure drop, to gravity and to interfacial 
shear. A Weber number We accounts for surface tension. With this rescaling, the velocity 
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profile for the undisturbed base flow in nondimensional form reads 
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( 3 ) 


0 < z < Sg- 


The ratio of pressure and gravity Reynolds number Rep/Reg results in a Fronde number, 
which, in combination with the pressure scale in Eq. ([^, represents a measure for the effect 
of applied pressure drop relative to gravity acting on the gas layer: 


Fr 


Vp 



( 4 ) 


To gain insight into the development of small disturbances ri{x,t) centred around the 
flat interface z = 0, we examine the linear stability of the interface in terms of a standard 
Orr-Sommerfeld-type analysis (a comprehensive formulation of this analysis can be found 
in Appendix ^ In this approach, the governing equations are linearized around the base 
flow, Eq. ([^, and inhnitesimally small perturbations in the associated streamfunction (j){z) 
are assumed to have wave-like solutions of the form ^{x,z,t) = (/{z), where a = 

ar + ictj is the complex wavenumber and u = + icu* is the complex angular frequency 

(both dimensionless understood). The imaginary part of these quantities, ctj and cuj, denote 
the spatial and temporal growth rates, respectively. If cjj > 0 the interface is considered 
temporally unstable and developing (initially) sinusoidal waves propagate with the phase 
velocity Vp = Ur/cHr- 

We solve the underlying generalized complex eigenvalue problem numerically by employ¬ 
ing a standard Chebyshev collocation method |1] (full description given in Appendix iQ, 
thereby adjusting the number of collocation points until convergence is reached. 


B. Nonlinear direct nnmerical simulation 


To capture and analyse the development of the flow beyond the linear regime, we use 
the two-phase flow solver presented in Reference |2I] for direct numerical simulation of the 
full Navier-Stokes equations. This in-house solver is level set method based and uses the 
continuum surface force (CSF) formulation to model surface tension effects [221123]. In level 
set formalism, the governing equations read as 


/ du 
— Vp 


u •VnI= 


Re, 


■V • [fi {Vu + Vn^)] + ( 0 ) niz, 


We 


(5a) 


V ■ u = 0 , 


d(t) 


+ u • 'V(j) = 


0 


(5b) 

(5c) 
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n = 




(5d) 


Here, 0 {x, t) is the level set function indicating the phase in which a point x lies (liquid 
phase for 0 < 0, gas phase for 0 > 0). Hence, the zero level set, (l){x,t) = 0, represents 
the interface r] (x, t). The level set function also determines the unit vector n normal to the 
interface and the interface curvature k in Eq. (5d). As the level set function differentiates 


between the two phases, it is also used to identify the respective density and viscosity through 
the expressions p = ifg (0) + r (1 — (0)) and /i = (0) + m (1 — (0))- The function 

(0) is a regularised Heaviside function, which is smooth in an interval [—e, e] around the 
interface with e set equal to 1.5 times the grid spacing. This interval also supports the 
regularised delta function <5^ (0) = (0) /d0. 

To discretize Eq. ([^, we use an isotropic MAC grid with a spacing that resolves the 
undisturbed liquid him with at least 30 grid points. Additionally, all simulations have 
been checked for convergence. On the implemented grid, vector quantities are dehned at 
cell faces and scalar quantities are dehned at the respective cell centres. A third-order 
Adams-Bashforth scheme is used to treat the convective derivative, while the momentum 
huxes are treated in a hux-conservative fashion employing a (semi-implicit) combination of 
Crank-Nicholson and third-order Adams-Bashforth method m- Pressure and the associ¬ 
ated incompressibility of the how are treated using standard projection method |25]. We 
use a combination of Jacobi’s method and successive over-relaxation on a red-black scheme 
to evaluate the predictor and corrector step. Moreover, a third-order (hfth-order accurate) 
WENO scheme [2S] is used to advect the level set function 0, which is subsequently reini¬ 
tialised applying a Hamilton-Jacobi equation and the algorithm formulated in Reference m- 
At the domain boundaries we apply standard no-slip and no-penetration conditions at the 
conhning walls, z = 0 and z = H (note that the coordinate system underlying the Eqs. ([^ 
is shifted by —5l compared to the one shown in Fig. [^, as well as periodicity in x-direction. 
The pressure is decomposed as p = p -|- (dP/dL) x, where p satishes the periodic boundary 
conditions in x-direction and dP/dL is a positive, constant, dimensionless pressure gradient. 


Solving the standard force balance in both phases (Eq. (A.l) and (A.5)) numerically gives 
the initial velocity held (base how). 

The solver is implemented in Fortran 90 using MPI to decompose the computational 
domain (Fig. in x-direction. This parallelization scheme makes efficient use of the ar¬ 
chitecture of the UK’s “Advanced Research Computing High End Resource” (ARCHER, 
http://www.archer.ac.uk) on which we run our high resolution simulations. 


III. TEMPORAL STABILITY ANALYSIS 

In this section we restrict ourselves to the study of the temporal evolution of an (initially) 
inhnitesimally small perturbation of the liquid interface. This temporal framework provides 
deep insight into the onset of interfacial waves and, thus, the complex dynamics of vertical 
hlms sheared by counter-current gas hows. It further enables us to map how regimes typical 
for such systems. 

We analyse the temporal stability for two distinct cases: 

• High density contrast: we assume a liquid density of = 1000 kg m“^, corresponding 
to a gas-liquid how typihed by an air and water combination. We demonstrate below 
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that the large density contrast leads to a complicated characterization of the instability, 
consisting of competing and coalescing linearly unstable modes. 

• Low density contrast: we assume a liquid density of pi = lOkgrn”^. The low-density- 
contrast scenario is popular in the simulation literature [25fl5U] for a number of rea¬ 
sons. In the present context, it corresponds to a system without mode competition. 
This provides a “clean” database of linear stability results, which can be used as an 
unambiguous benchmark for direct numerical simulations. 


For both considered cases, we further assume a liquid dynamic viscosity of pl = 500 ■ 10“® Pas 
and a surface tension of 7 = 1 ■ 10“^ N m“^, while the dynamic viscosity and density on the 
gas side are /xg = 10 ■ 10“®Pas and pc = lkgm“^, respectively. The flow is conhned by a 
channel of the height H = 0.01m. These values result in density ratios of r = 1000 and 10, 
a viscosity ratio of m = 50 as well as a gravity Reynolds number of Rcg = 313 and leaves 
the relative him thickness 5^ and the Fronde number Fr, which can be related to liquid and 
gas how rates, respectively, as the remaining parameters to determine the two-phase how. 

Throughout the linear stability analysis, i.e. both density ratio cases, we consider 5l & 
[0.02,0.14], whereas the Fronde number is varied within the interval Fr G [1.05,13] for the 
case of high density contrast. This corresponds to an absolute him thickness ranging from 
0.2 ■ 10“^ m to 1.4 • 10“^ m and an applied pressure drop dp/dx in the range of 10.8 Pa m“^ to 
1657.9Pam“^. Within this parameter space, we determine the temporal dispersion relation 
ujtemp = oj(^ar,C(i = 0) numerically on a grid with ASl = 0.005 and AFr = 0.25. The 
associated eigenvalue problem, Eq. (A.19), is solved for £ [0.05, ccc] with Aar = 0.05, 
where ac is the cut-oh wavenumber beyond which < 0. In the low-density-contrast 

case, we apply the interval Fr G [1.05,1.55], corresponding to dp/dx G [10.8, 23.6] Pam“^, 
with AFr = 0.025 and determine the dispersion relation for G [0.05,15] with Aar = 0.01. 
For each parameter set {d^^Fr) we further identify the pair that maximizes 

Ui as the linearly most unstable mode. 

For both density ratio cases, our analysis shows (Fig. that the temporal growth rate of 
that mode always attains positive values, which means the interface is inherently unstable 
within the investigated parameter spaces. However, the trends for the growth rate of the 
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FIG. 2. Temporal growth rate of the linearly most unstable mode in the parameter spaces 

of the two investigated cases, (a) r = 1000; (b) r = 10. 


most unstable mode are quite different in the two cases. Apart from a general increase of 











































(^temp increasing Fronde nnmber in the case with high density ratio, a valley stretching 
ont parallel to the horizontal axis at a Fronde nnmber of aronnd 2.25 becomes evident in 
the 5r-Fr-plane (Fig. |^). Moreover, a plateau-like region of nearly equal temporal growth 
rates appears at about Fr = 6.00 spanning all considered him thicknesses. Further insight 
into the nature of the system can be drawn from the length of the most unstable wave. 
Am = 27r/am’”^, and its extent relative to the respective him thickness, Xm/Si- Figure]^ 
shows the change of the latter quantity in the parameter space for density ratio r = 1000, 
revealing the complex interfacial dynamics of the how. Although the interface is susceptible 



FIG. 3. Wavelength Am of the linearly most unstable mode scaled by the corresponding relative 
film thickness dr in the parameter spaces of the two investigated cases, (a) r = 1000; (b) r = 10. 
The dashed line indicates zero phase velocity (onset of flooding) for reference. 

to a wide range of linearly most unstable waves, short-wave instability is predominant in 
this case, which is a consequence of the low surface tension considered herein [2]. 

The rich dynamics observed for the high-density-ratio case can be explained by looking at 
the dispersion relation of representative points in the parameter space (Fig. |^. It becomes 
apparent that multiple unstable modes are active, which are coalescing and competing with 
each other. It is this mode interaction that causes the complex structures presented above. 
In contrast, the case with low density ratio shows mostly one unstable mode throughout the 
considered parameter space. Even though further unstable modes are observed, these modes 
are very weak and appear only occasionally. Hence, for this case the system behaviour is 
less intricate (Fig. & and 1^). 

In order to identify the character of these unstable modes and hence the driving force of 
the instability in these scenarios, we characterize the energy transfer from the base flow by 
decomposing the disturbance kinetic energy into production and dissipation terms. However, 
as both presented parameter spaces are quite extensive, we restrict ourselves to the detailed 
analysis of the representative scenarios listed in Table |Tj The mode characteristics for all 
other parameter sets are deduced by inspection. 


A. Energy budget 

Following the approach in Reference [3], the rate of change of kinetic energy in both 
phases KINl,g can be decomposed into 

KINl + KINg = DISSl + DISSg + REYl + REYg + NOR + TAN, 


( 6 ) 
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FIG. 4. Selected dispersion relations for the case r = 1000. Solid lines: interfacial mode; dashed 
lines: shear mode in gas layer; dotted lines: shear mode in liquid layer; dot-dashed line: internal 
mode. Top to bottom: Fr G [2.5,3.0, 7.5, 8.5]; left to right: 6l G [0.03,0.08,0.13]. 


where all terms are spatially averaged contributions of the respective fluid to the overall 
energy budget of the disturbance. The terms DISSl^g represent energy losses due to viscous 
dissipation, whereas REYl g denotes wave-induced Reynolds stresses transferring energy 
between the base flow and the perturbation in the bulk of the two phases. Lastly, NOR and 
TAN describe the work done by normal and tangential stresses at the interface. The energy 
budgets for the investigated scenarios are given in Table |TT| where the individual terms have 
been scaled by the total rate of change of kinetic energy {KINl + KINq = !)• 

Decomposition of the energy budget in scenario TIC (Fig. |^) uncovers three unstable 
(active) modes: one interfacial mode as well as one shear mode in each fluid layer. The hrst 
mode (a = 40.05) originates from the viscosity and density contrast of the two fluids, which 
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r 

Scenario 

5l 

Fr 

Rep 

IFe 

temp 

Up,OS 

temp 

temp 

^i,DNS 


TIC 

0.13 

3.00 

939 

8.829 

40.05 

16.47 

1.6853 

- 

1000 

T2C 

0.08 

7.50 

2349 

55.181 

161.90 

0.82 

8.8447 

- 


T3C 

0.08 

3.00 

939 

8.829 

40.85 

5.97 

1.7662 

1.7783 


dT4C 


1.157 

362 

1.313 

3.99 

0.05 

0.3662 

0.3668 

10 

dT5L 

0.08 

1.179 

369 

1.363 

4.29 

0.00 

0.4669 

0.4634 

dT6F 

1.201 

376 

1.415 

4.59 

-0.05 

0.5829 

0.5790 


dT7F 


1.319 

413 

1.706 

6.19 

-0.30 

1.4347 

1.3806 


TABLE I. Conditions studied in detail using linear theory and nonlinear direct numerical simu¬ 
lations (case r = 10 only). The leading letter of the scenario name designates temporal stability 
analysis (prehx “d” for decreased density ratio), followed by a running number for this type of 
analysis; the trailing letter corresponds to the flow regime (Counter-current, Loading, Flooding; 
introduced in 


IIIB) in which the pair {6L,Fr) is situated. 


r 

Scenario 

a 

KINl 

KINg 

DISSl 

DISSg 

REYl 

REYg 

NOR 

TAN 


TIC 

1.30 

0.02 

0.98 

-0.01 

-1.39 

-0.04 

2.42 

0.00 

0.02 

1000 

TIC 

6.20 

1.00 

0.00 

-0.33 

-0.01 

1.31 

0.02 

0.00 

0.01 

TIC 

40.05 

0.51 

0.49 

-0.46 

-31.37 

-0.05 

-4.82 

-0.30 

37.98 


T2C 

69.50 

0.76 

0.24 

-0.31 

-18.65 

1.06 

-3.67 

-0.02 

22.57 


dT4C 

3.99 

0.06 

0.94 

-1.09 

-13.93 

0.00 

-1.77 

-0.46 

18.25 

10 

dT5L 

4.29 

0.09 

0.91 

-1.21 

-12.24 

0.00 

-1.19 

-0.44 

16.08 

dT6F 

4.59 

0.12 

0.88 

-1.36 

-10.67 

0.00 

-0.72 

-0.42 

14.17 


dT7F 

6.19 

0.32 

0.68 

-1.48 

-5.44 

0.00 

0.44 

-0.32 

7.79 


TABLE 11. Energy budgets of the active modes for the scenarios listed in Table |Ij 


give rise to velocity and stress disturbances tangential to the liquid interface driving the 
instability (“viscosity-gravity-induced instability”). Across the entire parameter space of the 
higli-density-ratio case, this interfacial mode is active (Fig. and, generally, its maximum 
growth rate as well as the associated wavenumber increase for increasing Fronde 
number yet slightly decrease with increasing relative him thickness, see Fig. and[^. 

The shear mode, or Tollmien-Schlichting mode, in the gas layer is also active for all 
values of above a minimum Fronde number ranging from 1.46 to 2.34 for increasing 
him thickness. Above this threshold, the mode becomes stronger with increasing gas how, 
though thicker hlms exhibit slightly lower growth rates. The opposite trend is observed for 
the wavenumber at maximum growth rate, which shifts towards lower values with increasing 
Fr and generally higher values for thicker hlms. Contrasting the behaviour of the shear mode 
in the gas layer, the liquid layer shear mode is only active for thick enough hlms {6l > 0.067). 
While this threshold value increases with increasing Fronde number, the maximum growth 
rate of this mode drops, thereby rendering it less unstable. Further complexity is added to 
the system dynamics by the fact that the energy contribution to this particular mode due to 
tangential stresses at the interface becomes larger for increasing Fr and reaches signihcant 
values {TAN k, 1). Eventually, TAN becomes the dominant term but with REY^ still being 
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of importance to overcome the restoring effects in both phase. This combination of energy 
sources driving the instability is also known as internal mode. 

Apart from this change in mode identity, there is a further internal mode emerging at 
Froude numbers ranging from 5.21 to 5.55 for increasing 6l, e.g. scenario T2C or Fig. |^. 
As this mode grows stronger with higher gas flow rates, the wavenumber at its maxi¬ 
mum growth rate shift towards lower values. Furthermore, the mode starts to coalesce with 
the dominant interfacial mode, which leads to the formation of a second “hump” in that 
dispersion relation (Fig. |^). For even higher Fr, the coalescing internal mode undergoes a 
change in identity itself and forms a second active interfacial mode (Fig. Ef and §)■ This 
“splitting” of the dispersion relation induces a region of rapidly changing wavelength of the 
linearly most unstable mode at intermediate Froude numbers (Fig. |^). Although mode 
coalescence occurs throughout the entire range of investigated him thicknesses, splitting of 
the dispersion relation is not observed for low values of 5^ (Fig. |^). It is further worth men¬ 
tioning that additional modes can become unstable, which are, however, mostly temporarily 
active and thus of minor signihcance. 

The multitude of active and coexisting modes leads not only to mode coalescence but 
also to a certain amount of mode competition. In the high-density-contrast case this mode 
competition mainly occurs at low Froude numbers, where the growth rate associated with 
the interfacial mode is still comparatively low. The shear modes at low and high values of 

on the other hand, exhibit stronger growth and therefore constitute the dominant mode 
(Fig. 1^ and |^). As Fr increases, the interfacial mode picks up strength and supersedes 
the shear modes as the dominating mode (Fig. and |^), which results in a jump of 
the wavenumber associated with the maximum growth rate. This jump explains the 

regions of relatively long waves at low and high values of 5^ in Fig. [^. 

As mentioned above, the low density-ratio-case exhibits mainly one unstable mode, which 
is due to the density and viscosity contrast of the two fluids (Table [n|). Although both dif¬ 
ferences account for energy transferred towards the disturbed flow, the contribution related 
to the viscosity contrast dominates. Hence, this mechanism is, in general, consistent with 
the so-called Yih mode [1]. It is further apparent that the relative fraction of kinetic en¬ 
ergy associated with the liquid phase increases with increasing Froude number. This rise, 
together with an enhanced energy dissipation, can be linked to more agitation in the liquid 
him as we will show in section |Vj The amount of energy dissipated in the gas phase, on the 
other hand, seems to drop, which is counter-intuitive for an increased Fr. Yet, dissipation 
in the gas does increase in absolute terms but at a slower rate as the total kinetic energy. 
That, in turn, leads to the seemingly decreasing rate of dissipation in the gas phase. The 
same effect can be seen for NOR and TAN. 

Another result worth noting is the change in sign of REYq as Fr increases, turning 
wave-induced Reynolds stresses from an additional dissipative energy “sink” to an energy 
“source”. This, in turn, can be explained by an unstable Tollmien-Schlichting mode ap¬ 
pearing in the gas stream (scenario dT7F) that delivers energy to the interfacial instability, 
thereby suggesting a transition to turbulence in the bulk of the gas phase. Therefore, 
this particular scenario can conceptually not be regarded as a “pure” Yih-type instability 
any more but tends towards an internal mode. This positive contribution of wave-induced 
Reynolds stresses to the instability occurs throughout the entire parameter space but is 
shifted to lower Froude numbers for thicker liquid dims. However, it has to be emphasized 
that the tangential stresses doing work at the interface are the dominant driving force of the 
instability in all presented scenarios for the low-density-ratio case (Table [H)). Across the en- 
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tire parameter space the maximum growth rate of this interfacial mode increases with 
increasing 5l and Fr (Fig. [^), whereas the associated wavenumber increases with 

increasing Fronde number but decreases for thicker dims, which in agreement with na. 
Overall, the system is predominantly receptive to long-wave instability under low-density- 
ratio conditions (Fig. [^). 


B. Flow regimes 

Fig. 0 shows the phase velocity Vp of the fastest growing wave developing on the interface 
for both density ratio cases. It becomes apparent that the parameter space is divided 
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FIG. 5. Phase velocity Vp of the linearly most unstable mode in the parameter spaces of the two 
investigated cases, (a) r = 1000; (b) r = 10. Zero phase velocity (dashed line) marks the onset of 
flooding. For comparison, the curve of zero interfacial velocity in the undisturbed base flow (dotted 
line). 

into two regimes: one in which developing waves exhibit a positive phase velocity and 
another in which the phase velocity is negative. With respect to the chosen coordinate 
system (Fig. [^, these regions correspond to waves propagating downwards and upwards, 
respectively. The vanishing Vp at the demarcation (dashed lines in Fig. between these 
regimes relates therefore to a standing wave. By crossing this line with increasing Fr, liquid 
adjacent to the interface begins to move upwards with the developing wave. This (partial) 
upward flow of the liquid phase is considered as flooding [31]. The curve Vp{6L,Fr) = 
0, which is referred to as loading curve in chemical engineering, indicates the operational 
limit of vertical counter-current gas-liquid flows. Thus, the region below the loading curve 
constitutes the counter-current regime. 

The flow map for the case with high density ratio (Fig. |^) reflects the complex interfa¬ 
cial dynamics described above. In regions with mode competition and changing dominant 
mode, the phase velocity changes at a great rate due to the jump of the wavenumber 
associated with the linearly most unstable mode. For low Fronde numbers, this leads to a 
region of flooding amidst the counter-current regime at the thin him end of the parameter 
space, while a region of comparatively slow waves occurs at the high him thickness end. 
In contrast, the phase velocity changes smoothly in the low-density-ratio case due to the 
absence of mode competition (Fig. [^). 














13 


To illustrate the effect of the linearly most unstable waves on the onset of hooding, Fig. 
also shows the curve of zero interfacial velocity (dotted line) for the undisturbed base flow 
of the respective case, which can be seen as the analogue to the actual loading curve. In 
the high-density-ratio case the two curves (with exception of the “island” region) almost 
coincide, i.e. linear waves have no significant influence on the loading curve (Fig. 

This can be explained by looking at the expression for the phase velocity, which can be 
decomposed (trivially) as Vp = uo^int + ui{a), where Mo,mt denotes the interfacial velocity 
of the undisturbed base flow. Yet, a good approximation for ui{a) can be found in the 
dispersion relation for free-surface capillary waves [nii32], such that ui{a) ~ Ucap/c(, where 
^cap — / {pL + Pg) (dimensional understood). Even in the absence of this particular 

approximation, it can be argued that Ui{a) should be symmetric in the densities (at least 
for the present vertical configuration), and should be a monotonically decreasing function 
of the densities, as generally the inertia tends to affect the phase speed in a reciprocal 
manner, i.e. the larger the inertia the smaller the phase speed. Thus, it can be seen that 
that Vp —)■ Uo^int for high density ratio. This also implies that interfacial waves do alter the 
hooding curve in the case with low density ratio. The counter-current regime is significantly 
reduced by the instability for large parts of the parameter space (Fig. |^). However, below 
a threshold of 6l = 0.036 the onset of hooding is pushed towards higher Fronde numbers. 


IV. ABSOLUTE AND CONVECTIVE INSTABILITY OF THE LINEAR 

DYNAMICS 


Complementing the temporal analysis of the interfacial instability, we study the response 
of the system to a pulse on the liquid interface. In particular, we are interested whether this 
initially localized disturbance spreads both up- and downstream and eventually contaminates 
the entire how, rendering the system absolutely unstable. In case the disturbance propagates 
away from its source, the how is considered convectively unstable. 

To determine the nature of the instability in this context, the complex dispersion relation 
D (a,a;) = 0, which is obtained by solving the eigenvalue problem of Eq. (A. 17) numerically 
for a range of complex wavenumbers a = ar + idi, has to be evaluated against conditions 
essential for absolute instability as outlined in Reference [7j: (i) a positive imaginary part of 
the angular frequency := Ui (as) > 0 at a saddle point as in the complex a-plane, where 
as solves dcu/da = 0, forms the necessary condition; (ii) to satisfy the sufficient condition, 
spatial branches a^ {u) that originate from opposite halves of the a-plane have to coalesce 
at the saddle point as, forming a genuine pinch point (this coalescence corresponds to 
the formation of a branch point/cusp at in the complex cj-plane |33]). Meeting both 
conditions will result in growth of the disturbance at its source with the absolute growth rate 


Although the described procedure seems straightforward, inspection of the results of 
a spatio-temporal Orr-Sommerfeld (ST-OS) analysis requires great care in order to avoid 
misinterpretation due to the complexity of the dispersion relation that can arise from the 
multivalued nature of the eigenvalue problem as well as specifics of the investigated problem, 
such as applied boundary conditions or multiple unstable temporal modes (see Fig. |^. 
Hence, to conclude the character of the instability correctly, it is necessary to examine the 
global topography of frequency over wavenumber for each relevant set of parameters. As we 
are interested in identifying convective/absolute instability (C/A) transition, the described 
method is not practical given the large parameter space considered herein. Instead, we make 
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FIG. 6. Global topography of Wj for scenario ST4 (r = 10, 6 l = 0.08, Fr = 1.201, Rcp = 376, 
We = 1.415) as obtained from spatio-temporal Orr-Sommerfeld analysis. The cross marks the 
pinching saddle point. 


use of the analytical connection between temporal and spatio-temporal frequency presented 
in Reference [34] . 

The authors derived an expression for Ui (a^, w) based on analytical continuation of 
purely temporal stability properties, specihed on the real axis, into the complex a-plane. 
Accounting for up to second-order terms, the quadratic approximation (QA) of the imaginary 
part of the complex frequency reads 




uJi {ar,ai) = (ar) + 


dar 


-a,; 


da? 


rj 2 

a,. 


(7) 


Similarly, we approximate (a^) in the vicinity of the temporally most unstable mode by 
second-order polynomials. With the conditions dui/dar = duijdai = 0 applied to Eq. ([^, 
we can readily estimate the position of the saddle point as and the absolute growth rate 
for a given set of flow parameters, requiring results from temporal linear stability analysis 
only. Even though this approach seems to circumvent the pitfalls associated with the full 
spatio-temporal linear stability analysis outlined above, one has to keep in mind that the 
quadratic approximation is, strictly speaking, only valid inside a disc of convergence with 
radius R, where R is the minimum distance from the centre of the complex Taylor series 
{ar = cd^P, 0) to the nearest singularity of uj (a). 

Like in the purely temporal framework, we study several scenarios in more detail (Ta¬ 
ble III). Using spatio-temporal Orr-Sommerfeld analysis and quadratic approximation as 
outlined above, we determine the saddle point as as well as the corresponding absolute 
growth rate uji^ for all listed scenarios and further perform direct numerical simulations for 
the scenarios STl and ST4. In the following, we want to compare and discuss the obtained 
results for scenario ST4 in more detail. 

As shown earlier already, the global topography of oji in the complex a-plane (Fig. 
is rather complex. First, the dispersion relation contains two saddle points, of which both 
may be dynamically relevant. The conhnement of the flow by the channel walls has, fur¬ 
thermore, created a discrete pole on the imaginary axis [Mj, {ar,ai) = (0,3.34), which has 
implications on the character of the saddle point closer to that particular singularity. Lastly, 
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Scenario 

Sl 

Fr 

Rcp 

We 

Method 

Oir,S 

Cti,S 

<Xi,0 

STl 

0.08 

1.075 

337 

1.134 

QA 

DNS 

3.34 

-2.24 

-0.0877 

<0 

ST2 

0.08 

1.157 

362 

1.313 

QA 

ST-OS 

4.01 

4.02 

0.24 

0.24 

0.3632 

0.3632 

STS 

0.08 

1.179 

369 

1.363 

QA 

ST-OS 

4.35 

4.43 

0.66 

0.65 

0.4422 

0.4423 






QA 

4.72 

1.02 

0.5205 

ST4 

0.08 

1.201 

376 

1.415 

ST-OS 

4.87 

0.97 

0.5221 






DNS 

- 

- 

0.5269 


TABLE III. Scenarios with low density ratio (r = 10) studied in the spatio-temporal (ST) frame¬ 
work. 


the multivalued nature of the dispersion relation becomes apparent by the branch cut just 
below the real axis. Although these features make the hnal characterisation of the insta¬ 
bility more difficult, the saddle point at {ar,ai) = (4.87,0.97) clearly appears as a result 
of the coalescence of spatial branches emanating from opposite half-planes. It is therefore 
also pinch point and contributes to spatio-temporal growth at a rate of ujifi = 0.5221. The 
positive value of o;* indicates that the disturbance travels upwards, which is conhrmed by 
DNS (Fig. [^). On the other hand, the saddle point at {ar, ai) = (0.62,0.65) is not a pinch 
point. In fact, the spatial branch (cu) (above the saddle point) is a closed curve, whereas 
the a~ {uj) branch does not (entirely) originate from the negative half-plane. Hence, this 
saddle point is dynamically irrelevant regarding absolute instability. 

As mentioned before, the singularity closest to the position of the temporally most un¬ 
stable mode {ar = 0 ) determines the disc of convergence in which the quadratic ap¬ 

proximation can be applied with conhdence. In scenario ST4, the conhnement pole at 
{ar, Qfj) = (0, 3.34) limits the outermost radius R of this disc to about 5.68. Equation ([^ is 
thus convergent across the relevant section of the complex a-plane depicted in Fig. The 
approximated position of the pinching saddle point as well as the corresponding growth rate 


agree well with spatio-temporal OS analysis, see Table III 


We further compare the impulse response for the parameters of scenario ST4 captured by 
DNS against linear theory (ST-OS). In contrast to simulations within the purely temporal 
framework, a Gaussian pulse with a height of 1 ■ 10“^ and a standard deviation of 1 ■ 10“^ 
is implemented on the interface as an initial condition to study the spatio-temporal nature 
of the flow. Growth of this perturbation at its source then constitutes absolute instability. 
However, the developing disturbances will inevitably contaminate the pulse source due to 
the implemented streamwise periodic boundary conditions. To delay this contamination 
sufficiently, the channel length is set to = 20. Using the pulse norm 


n {x, t) = 


\w {x, z, t)f dz 


1/2 


the space-time plot in Fig. shows the temporal evolution of the perturbations along the 
streamwise coordinate. It can be seen that, starting from its initial position at x = 16, 
the pulse is convected upwards, which is in accordance with the results from linear theory 











16 


' 


LJ4.1E-11 

0 2.5 5.0 7.5 10.0 12.5 15.0 



FIG. 7. (a) Space-time plot of the norm n {x, t) for scenario ST4; (b) Comparison between growth 
rate at the source of the disturbance obtained from DNS and spatio-temporal Orr-Sommerfeld 
analysis. 


mentioned above. This upward motion triggers a pressure disturbance moving ahead of 
the pulse. However, this “shock” decays sufficiently before it re-enters the computational 
domain, thus avoiding early contamination of the instability source (see left-hand side of 
Fig. 1^). Further information from this plot is extracted in Fig. Izf. where the growth of 
the instability at its source is given. After an initial transient period, the growth rate of 
the instability stabilises and follows the value predicted by ST-OS and QA closely, therefore 
validating the DNS results. Thereafter, a second pressure shock develops {t ~ 11.5) and 
travels upwards to eventually contaminate the pulse source at around t = 14.0. 

The excellent agreement that has been established between linear theory (ST-OS) and 
quadratic approximation in ST4 can also be observed for the scenarios ST2 and STS (Ta¬ 
ble III). For scenario STl we were not able to determine the absolute growth rate of the 
system by means of OS analysis and saddle point method due to the multivalued nature of 
the associated eigenvalue problem. In fact, the dynamically relevant saddle point appears in 
the lower half-plane below the branch cut arising near the real axis. To analyse this saddle 
point a laborious reconstruction of the corresponding part of the Riemann surface from dis¬ 
crete eigenvalues would have to be carried out. However, we avoid this procedure by using 
the quadratic approximation, which indicates scenario STl to be convectively unstable. A 
further direct numerical simulation conhrms this (not shown). 

Given these results, we apply the QA to identify the spatio-temporal nature of the sys¬ 
tem throughout the wide parameter space considered in the low-density-ratio case. The 
calculated growth rates Ui^ are shown in Fig. where the dashed lines demarcate the C/A 
boundaries. It becomes apparent that the system is absolutely unstable in almost the entire 
domain with exception of two narrow bands at the low Froude number and low him thickness 
end of the parameter range, respectively. In view of these results and those of § [III B[ it 


appears that C/A transition and the onset of hooding are not closely related for the case of 
density ratio r = 10. 

Regarding the accuracy of these results, it has to be mentioned that the relative er¬ 
ror (cuqa — tJsT-os) / |i^ST-os| between QA and linear theory increases with decreasing him 
thickness and increasing Froude number. At the same time, the quadratic approximation 
generally underestimates the “true” value of Ui,o as obtained by ST-OS, which leads to a 
larger regime of absolute instability (towards thin hlms) than displayed in the above hgure. 
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FIG. 8. Absolute growth rate Uifi at the saddle point (contours) as obtained by analytic continua¬ 
tion (quadratic approximation) from purely temporal quantities. The dashed lines demarcate the 
transition between convective and absolute instability. 


Finally, the mode coalescence and mode competition encountered in the case with high 
density ratio (e.g. Figure makes it difficult to identify the nature of the instability in 
a spatio-temporal framework with both the present methods. To study this manifestation 
of the instability, alternative techniques, such as linearized DNS [ 6 ] (wherein the linearized 
equations of motion are solved as a Cauchy problem, but still within the framework the 
Orr-Sommerfeld linear operators), might be more suitable. A detailed investigation of this 
particular regime, i.e. high density contrast, will be left to future work. 


V. NONLINEAR WAVE DYNAMICS 


As temporal linear stability analysis permits the analysis of inhnitesimally small interface 
perturbations only, we carry out direct numerical simulations to study the evolution of these 
disturbances up to hnite amplitudes for both low and high density contrasts. To that end, 
we use the in-house solver described in 


and initial interface location 


IIB with streamwise-periodic boundary conditions 


N 

7] {x,t = 0) = 6l + cos {anX + (fn) , ( 8 ) 

n=l 


where Aq is some small initial amplitude (herein Aq = 1 ■ 10“^ ), N is the number of linearly 
unstable modes initialized, an = n {271/L^) n is a wavenumber in streamwise direction and (fn 
is a random phase. Even though periodic boundary conditions do not reflect the behaviour 
of a real system, it is appropriate to consider this setup as it allows for a rigorous comparison 
of DNS results with linear theory as well as unambiguous results of the Fourier transform 
taken of the interfacial wave at hnite times | 21 j . 
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A. Low density contrast 

To allow for rigorous comparison with linear theory, we perform direct numerical simula¬ 
tions for each flow regime of the low-density-contrast case (lower part of Table [^. For that 
purpose, we set the wavenumber in Eq. (|^ equal to the linearly most unstable mode 
of the respective scenario. Figure shows the L^-norm of the wall normal velocity 
perturbation w over time for scenario dT4C and it can be seen that the disturbance grows 
with the rate predicted by Orr-Sommerfeld theory. Comparing the streamfunction 




FIG. 9. Comparison of nonlinear DNS with linear theory through semi-analytical Orr-Sommerfeld 
analysis for scenario dT4C. (a) shows the wave growth rate calculated as the L^-norm of the 
perturbation rc-velocity; (b) is a plot the streamfunction normalized by the wave height at the 
position of the wave crest at t = 4.2. 

(pj (z) at the crest of the developing wave also yields excellent agreement (Fig. [^). Also 
the other scenarios follow the theoretical predictions in terms of growth rate in an equally 
convincing fashion (Table [T]). 

It can further be seen in Fig. that the wave initially grows exponentially before non¬ 
linear effects gain importance at about t = 8.0. Beyond that point, wave growth slows down 
and eventually saturates, leading to a steady nonlinear wave of constant amplitude travel¬ 
ling along the interface. Figure [TO^ depicts the early stage evolution of the interface up to 
saturation for the counter-current scenario. It can be appreciated that the wave developing 
on the interface moves downwards as is predicted by linear theory. More insight into the flow 
features of the counter-current flow scenario is given by Fig. [TT^ . Plotted in a hxed frame of 
reference, a large anticlockwise rotating recirculation zone, positioned on the “downwind” 
side of the wave, is revealed in the gas phase. This vortex is sandwiched between the main 
flow of the gas and a thin layer of gas along the entire length of the interface, which is 
dragged downwards by the liquid due to interfacial shear stress. The liquid phase itself does 
not exhibit any major disturbances. 

According to linear theory, the loading point of the system is reached by increasing the 
Froude number to Fr = 1.179 (scenario dT5L, see Fig. [^). At this point, the phase velocity 
of the wave vanishes as gravitational and lift forces on the wave hump balance each other, 
resulting in a standing wave. The evolution of this wave in its early stage is depicted in 
Fig. [T0|3 . It can be apprehended that the initial upward motion of the wave slows down after 
it is saturated and eventually comes to a complete halt around the position last indicated in 
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FIG. 10. Typical evolution of the disturbed interface in the three flow regimes (low-density-contrast 
case) within the time interval t = [0,17.5] at steps of At = 2.5 (dotted lines). The solid lines show 
the inital and final shape of the interface, respectively, (a) counter-current flow, dT4C; (b) loading, 
dT5L; (c) flooding, dT6F. 



FIG. 11. Pressure perturbation field, streamlines and liquid interface in a wall-fixed frame of 
reference for the scenarios near the loading curve (low-density-contrast case), (a) counter-current 
flow, dT4G, t = 27.7; (b) loading, dT5L, t = 34.6; (c) flooding, dT6F, t = 23.3. 


the plot. As the system evolves further the saturated wave begins to move downwards with 
the bulk of the liquid (not shown here). This is attributed to nonlinear effects which lead 
to an anticipated deviation from the system behaviour predicted by linear theory. However, 
it is worth mentioning that the deviation is relatively small [vp = 0.0103) and, hence, linear 
theory gives a good indication of the development and behaviour of the liquid interface for 
the loading scenario (dT5L). Moreover, the increased Fronde number leads to a decrease of 
the mean liquid streamwise velocity, which, in turn, reduces the extent of the thin layer of 
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downward gas flow along the liqnid interface. Thus, the gas side vortex developing in the 
wave trough moves closer to the interface. The enhanced flow velocity on the gas side, on 
the other hand, stretches this eddy along the entire trough region (Fig. [T^). On the liquid 
side, contrary to the counter-current scenario, a small anticlockwise rotating vortex occurs 
in the wave body, which indicates a negative (upwards) streamwise velocity of the liquid at 
the interface in the vicinity of the wave crest. However, the bulk of the liquid, which refers 
to the region of the film between channel wall and wave trough, remains largely unaffected. 

By increasing the Fronde number beyond the loading point, as in scenario dT6F, lifting 
forces on the wave body due to the enhanced gas flow overcome the gravitational forces, 
causing interfacial waves to travel upwards (Fig. [T^). Similar to the previous scenarios, 
the developing wave saturates and forms a coherent structure. However, with higher Fronde 
numbers the general trend of increasing growth rate and decreasing amplitude becomes 
apparent. Furthermore, the change of direction in which the interfacial wave propagates 
leads to a more agitated liquid film, especially in the wave body (Fig. 11:). Compared 


to scenario dT5L, an extended region of liquid near the wave crest experiences upward 
movement with the wave body and becomes recirculated in an enlarged eddy. The bulk of 
the liquid, on the other hand, essentially remains unaltered also in this scenario. 

After having discussed the influence of the Fronde number, i.e. gas flow rate, on the 
gas-liquid flow in a qualitative manner, we now want to turn our attention to a qualitative 
description of the nonlinear wave dynamics. To understand the genesis of these, we compute 
the spectra of the interface height for each time step. The spectra of the first three harmonics 
are shown in Fig. 12 r. Matching the rate given by Orr-Sommerfeld theory, the first harmonic 
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FIG. 12. Comparison of direct numerical simulations with weakly nonlinear theory for the counter- 
current scenario (dT4C). (a) interfacial spectra; (b) Stuart-Landau equation, eq. , fitted to 
fundamental mode. 

(«! = = 3.99) grows exponentially fast at the beginning. At the same time, the higher 

harmonics (02 = 7.98,03 = 11-97) are linearly stable, as predicted by the same theory. 
However, these harmonics eventually also undergo exponential growth, whereat the 
harmonic grows at a rate («o)- As time goes by, the growth rate of all harmonics 

decreases and the amplitudes saturate simultaneously. It is thus apparent that the dynamics 
of the higher harmonics are “slaved” to the first harmonic. This temporal development of the 
interfacial spectra is perfectly consistent with weakly nonlinear theory HI. The framework 
of the same theory allows to model the temporal development of the first harmonic using 
the Stuart-Landau (SL) equation |35l|36]. Including up to quintic-order terms, the absolute 
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value of the finite amplitude follows 


d 

dt 


l^«i 


hi 1^ ai f + h2 1^ ai 1^ + hs 1^ ai I : 


(9) 


where pi is twice the temporal growth rate and ^ 2,3 is twice the real part of the 

Landau coefficients |H|. These coefficients have been fitted numerically for the different 
scenarios discussed above and are listed in Table IV together with the root mean square 
deviation (RMSD) and R^-value of the best fit, illustrating the excellent agreement between 
theory and direct numerical simulations. The negative value of p2 confirms the existence 


Scenario 

hi 

h2 

h3 

RMSD 

R2 

dT4C 

0.7344 

-867.1 

252355 

1.1306 • 10-3 

0.9932 

dT5L 

0.9290 

-913.9 

0 

8.1751 • 10-^ 

0.9933 

dT6F 

1.1582 

-1410.1 

0 

5.7185 • 10-^ 

0.9957 

dT7F 

2.7613 

-5466.6 

0 

9.0486 • 10“^ 

0.9749 


TABLE IV. Landau coefficients (real part, fitted to DNS results) for the scenarios investigated 
under low density contrast. 

of the supercritical bifurcation observed in all scenarios. For scenario dT4C, the relatively 
large value of ps points towards a longer transition phase with decreasing growth of the 
wave amplitude leading from the initial stage of exponential growth to saturation in this 
scenario. In contrast, the wave in case dT5L, dT6F and dT7F reaches its saturated state 
faster, which is reflected by the vanishing of ps in Eq. (§. The comparatively low R^-value 
in scenario T7F can be explained by an overshoot of the amplitude in the transition phase 
due to transient effects before the wave settles at a stable equilibrium amplitude. 


B. High density contrast 


To gain insight in the nonlinear wave dynamics under high-density-contrast conditions, 
we perform direct numerical simulations with the system parameters corresponding to 
Fig. 1^, scenario T3C (Table [^. This system exhibits two distinct linearly unstable modes, 
one long-wave shear mode and a short-wave interfacial mode, which are accounted for by 
dREYc = 1-549 and ax an = 40.27 in Eq. ([^, respectively. It is due to the short-wave na¬ 
ture of the interfacial mode combined with the high inertia of the liquid phase that the flow 


characteristics near the interface (Fig. 13) differ significantly from those observed in the low- 
density-contrast case. Unlike the scenarios presented in the previous section, no prominent 
features, such as recirculation zones, emerge within the wave train. However, a vortex layer 
with the familiar “cat’s eye structure” develops at the demarcation between the bulk of the 
gas and a thin gas layer dragged downwards by the liquid due to interfacial shear stress. 
Thereby, the vortices are pinched between two consecutive high-pressure regions, which are 
forming on the downstream side of the short-wave crests. Further snapshots indicate that 
this vortex layer is unstable to secondary instability. Hence, the wave form in Figure 13 
should not be regarded as a quasi-steady state. 

Similar to the low-density-ratio case, the nonlinear dynamics of the interfacial mode 
appear to be consistent with Stuart-Landau theory, albeit these dynamics develop faster 
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FIG. 13. Pressure perturbation field, streamlines and liquid interface in a wall-fixed frame of 
reference for the scenarios near the loading curve (high density contrast), (a) t = 1.08; (b) t = 1.71; 
(c) t = 2.025. 
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due to the higher growth rate of the hrst harmonic {ax an = 40.27). Evidence is provided 
in Figure 14 Consistency with the Stuart-Landau theory is especially clear in the time- 




FIG. 14. Spectral analyses confirming the relevance of the Stuart-Landau theory for the high- 
density-contrast case, (a) Spectra in the spatial domain showing the evolution of relevant wavenum¬ 
ber normal modes; (b) Power-spectral density /(w) in the frequency domain. 


frequency domain. In particular, in Figure [14)3, a power-spectral density 


/ M = 


w (xo, Sl, t) e dt 


xo = 0.007 


( 10 ) 


is shown (T corresponds to the duration of the simulation). A well-dehned global maximum 
is observed at at Ur = 239.8, corresponding to the frequency of the most-dangerous TAN 
mode at a = 40.27 in the spatial domain. Successive maxima at cUr = n(239.8) with 
n = 2,3, ■ ■ ■ indicate that the “overtones” of the TAN mode are slaved to the TAN mode 
itself. 

The maxima in the power-spectral density function, although well dehned, are by no 
means sharp. The broadening of the maxima here is a sign that the wave train is not strictly 
periodic, and that a very large number frequencies is present in the dynamics. Crucially, the 
broadness of the lines is a function of the simulation time: for longer simulations, more and 
more frequencies come into play (albeit that the same maxima remain dominant throughout). 
This indicates the onset of chaos. Thus, the Stuart-Landau theory manifests itself as the 
leading-order approximation to a chaotic dynamics, albeit that the wave train in Figure 1^ 
is subject to secondary instability. 

Other second-order effects are key to understanding the secondary instability. The vortex 

In particular, a long- 


layer in Fig. 13 r is not steady but breaks down at later times (Fig.[l3[: 
wave perturbation to the vortex layer (wavelength on the domain scale) is seen to coincide 
with the signihcant growth of the long-wave shear mode. Therefore, it would appear that 
a complicated secondary instability sets in, involving a destabilization of the vortex layer 
by perturbations that are fed by the long-wavelength linearly unstable (but subdominant) 
mode. At this stage, individual waves steepen even further, up to the point where wave 


overturning is imminent (Fig. 13 d 


One may furthermore recall that the strong presence of a shear mode in the system is 
a signal of the supercritical transition to bulk turbulence in one of the phases. This is not 
directly relevant in the present scenario, as the secondary instability of the laminar wave 
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train sets in before this transition. In either case, it is clear that secondary instability may 
inhibit the operation of the system in a qnasi-steady laminar state at high density ratios. 


VI. CONCLUSIONS 

We have presented a comprehensive stndy on two-dimensional laminar flow of a vertical 
liqnid him sheared by laminar connter-cnrrent gas in a conhned channel. This stndy tries 
to fnrther elucidate the nature of interfacial instability in such two-phase hows using several 
complementary methods, namely Orr-Sommerfeld analysis, energy budget analysis as well 
as high resolution direct numerical simulations. Two sample systems have been selected for 
investigation: one with a high density contrast {pi/pc = 1000) and a second with a low 
density contrast {pl/pg = 10). In both cases, the same viscosity contrast {pl/pg = 50) 
and comparatively low surface tension (7 = 1 - 10“^Nm“^) was used. In our study we 
focussed on analysing the inhuence of liquid him thickness and applied pressure drop on the 
development of interfacial waves. 

Temporal linear stability analysis reveals that the liquid interface is inherently unstable 
for both cases. In the system with high density ratio short-wave instability is predominant, 
whereas the low-density-contrast case tends to favour long-wave instability. Furthermore, 
the instability is governed by a multitude of coexisting unstable modes (interfacial mode, 
shear mode in both phases, internal mode) under high-density-ratio conditions, where the 
latter two modes indicate the onset of turbulence in the bulk of one of the phases. In 
contrast, the low-density-ratio system exhibits only one unstable mode, which is consistent 
with the Yih mode [1] and a manifestation of the viscosity contrast. 

Additionally, we use the phase velocity of the linearly most unstable mode to identify two 
distinct flow patterns: counter-current flow and flooding. In the former regime, the phase 
velocity is positive, which relates to interfacial waves propagating downwards. By contrast, 
flooding is characterised by upward travelling waves, leading to partial flow reversal in the 
liquid him. During the transition from counter-current to hooding regime, standing waves 
evolve at the loading point. In the case with high density contrast, we found a region of 
hooding amidst the counter-current regime, which is the result of mode competition (shear 
mode in the gas layer vs. interfacial mode). 

We have further determined the nature of the instability in the spatio-temporal framework 
under low-density-ratio conditions. Using this information, we established a second how map 
indicating the transition from convective to absolute instability (C/A). Besides standard 
Orr-Sommerfeld analysis, we also adopted the analytical connection between temporal and 
spatio-temporal growth rates in the linear regime as presented in Reference [M]. This 
approach, which is based on analytical continuation, circumvents possible difficulties in 
identifying the absolute growth rate that are associated with the multivalued nature of the 
eigenvalue problem and specihcs of the problem at hand. Compared with OS analysis and 
DNS, this method shows good agreement and is therefore appropriate to accurately estimate 
the absolute growth rate of the instability. We hnd that the system is absolutely unstable 
in most parts of the parameter space considered herein with the exception of two narrow 
bands at low applied pressure drop and low him thickness, respectively. In the high-density- 
ratio case, mode competition and mode coalescence between the multiple unstable modes 
hindered the mapping of convective and absolute instability in the respective parameter 
space. For such systems, linearized DNS [6] may allow for more conclusive results. 

To assess the development of interfacial waves up to hnite amplitudes, we perform direct 
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numerical simulations for parameters within the established flow regimes of the low-density- 
ratio case, using a level set method based solver that has been developed in-house. These 
simulations show excellent agreement with linear theory during the stage of exponential wave 
growth and conhrm the determined flow patterns. The simulations also show saturation 
of the waves once nonlinear effects become important. A Fourier analysis reveals that the 
growth of the higher harmonics of the interfacial waves is coupled to that of the fundamental 
in a fashion which is consistent with weakly nonlinear theory. The growth of the hrst 
harmonic also agrees well with the Stuart-Landau model, thus underpinning the weakly 
nonlinear nature of the instability and, moreover, suggesting the existence of a supercritical 
bifurcation. Regarding high-density-contrast conditions, the dynamics of the (short-wave) 
interfacial mode appear to be similar to those observed under low density ratio. However, 
direct numerical simulations suggest that the emergence of an additional (long-wave) shear 
mode triggers a secondary instability, which leads to a chaotic wave train showing signs of 
imminent wave overturning. 

In summary, the combination of generic complementary (semi-)analytical and numerical 
methods presented herein yields a comprehensive and convincing characterization of inter¬ 
facial instability in vertical counter-current gas-liquid flows. We are therefore conhdent that 
this rigorous methodology can be employed to further elucidate the dynamics of parallel 
shear flows in a wide range of technically relevant parameter regimes, such as flows with 
high viscosity contrast or non-steady gas flow. Beyond that, the outlined approach may be 
used as a starting point for the future study of heat and mass transfer phenomena in vertical 
gas-liquid flows. 
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Appendix: Pull formulation of the linear stability analysis and numerical solution 

1. Base flow velocity profile 


Here, we give a detailed formulation of the governing equations underlying the linear 
stability analysis undertaken in §§ HI and IV as well as a description of the numerical 


method used to solve the corresponding generalized complex eigenvalue problem. 

Under the assumption of a steady, spatially uniform, laminar and incompressible flow 
in both phases, the Navier-Stokes equations governing the undisturbed base flow reduce 
to standard balances between pressure as well as viscous and gravitational forces. For the 
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liquid film, this balance can be written as 

d^ho dp 

hL-Ti— n—^ PlQ — 0, 
az^ ax 


-di < z < 0, 


(A.l) 


where Uq is the dimensional base flow velocity in the respective phase (tildes denote di¬ 
mensional quantities). Equation (A.l) is subject to no-slip at the liquid side channel wall, 

(A.2) 


2 = —hi, and continuity of tangential stress at the interface, z = 0: 

duo 


uo {z = -di) = 0, pl 


d2 




2=0- 


Thus, integration of Eq. (A.l) yields 

“0 {z) = ^ f - PL^ {z^ -dl) - — {z + di), -di <z<0 (A.3) 

2fiL \dx / ' Pl 

as the velocity prohle for the liquid him. The interfacial velocity, which constitutes one of 
the boundary conditions of the gas layer, reduces to 


^0,int ^0 (^ 0 ) 


1 


dp ^ \ j2 Tni 




(A.4) 


The velocity prohle for the laminar gas layer is derived analogous to the liquid layer. We 
therefore write the force balance as 

- ^ + PgP = 0, 0 <z <dG, (A.5) 

which is subject to continuity of velocity and shear stress at the interface: 

dho 


^0 (z 0) Uo^inti pG' 


d^ 




2 = 0+ 


Applying these interfacial condition on Eq. (|A.5|) yields 

1 


^0 \Z) T 


dp \ _2 ~ n ^ ~ ^ J 

o V h- PGg]z -2;, 0 < ^ < Pg 

2pG \dx J pg 


(A.6) 


(A.7) 


as the gas side velocity prohle. To determine the interfacial shear Tint we apply the no-slip 
condition at the gas side channel wall, Uq {z = dc) = 0: 


dp 


'^int 


T - pL9]dL - dL + 7^— -7-PgP Mg-M = 0. 


dp 


'^int 


2pL Vdx J Pl 2,PG \dx / “ Pg 
In summary, the velocity prohle of the undisturbed base how in dimensional form reads 

^ 1 


(A.8) 


2^, VS-'’"») (5^ - 4) - (i + . 


—di < z <0, 


“0(^) = < - 


dp \ 2 Ant , 

2^, 


(A.9) 


, 1 / dp \ 2 Ant ~ n ^ ~ ^ ^ 

+ 7^— n- Pig ] z - 2 ;, 0 < 2 ; < dG. 

2pG \dx J Pg 


Applying the nondimensionalization scheme of Eq.([^ and ([^ hnally leads to the dimen¬ 
sionless form of the base how velocity prohle as presented in Eq. ([^. 





















27 


As mentioned in 


2. Perturbation equations 

TTM we introduce a small disturbance that shifts the flat interface 


from z = 0 to z = T], where |? 7 | ^ 1. This (dimensionless) wave elevation gives rise to 
perturbations in the flow held of the form: 


(u, w, p) = {uq (z) + 6u (x, t ), 6w (x, z, t ), po (z) + 6p (x, z,t)), 


(A.10) 


where the subscript zero denotes base how quantities and the 6 quantities are inhnitesimally 
small perturbations. We use these variables of the how held and obtain the linearized 
Navier-Stokes equations, which, after elimination of the pressure, further yield the linearized 
equation for the wall-normal vorticity Suz in both phases {j = L,G): 


d 


d 


rrii 


dt dt Vi Re, 


-V^ ] 6uz = -6w-^Uo, 


(A.ll) 


with {rL^rc) = 1), {mi, me) = {m,l) and 


6uy = -^6u — -^6w. 

oz ox 


(A.12) 


It is further convenient to use the streamfunction representation (5 m, 5w) = {d^/dz, —d^/dx) 
for the two-dimensional disturbance velocity held, hence 


5uz = V"<I>. 


(A.13) 


Assuming a wave-like solution of the form $(a;, ^,f) = e^^“^ the streamfunction. 

Equations (A.12) and (A.13) lead to the Orr-Sommerfeld equations governing the stability 


of the liquid interface: 


la 


Uq 




a 


d?Uo 

dz^ 


m,- 


/d2 


- 




a 




(A.14) 


where a = Or + ia* and ca = + itu* are the complex wavenumber and angular frequency, 

respectively. This equation is subject to no-penetration and no-slip conditions at both 
channel walls: 

0 {-h) = {-Sl) = 0 (5g) = {Sg) = 0. (A.15) 

dz dz 

Furthermore, conditions for continuity of velocity and tangential stress as well as a jump in 
the normal stress due to surface tension are applied to match the streamfunction across the 
interface at z = 0 (we use the notation c = oj/a): 


<pL = (j)G, 


(A.16a) 


d(pL d(pG 
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/ d^ 2 






0L 

( dMo 

dMo 

C - Mo 

^ dz 

0+ dz 
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f d^Mo 

C — Mo V dz^ 
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dz^ 


(A.16b) 
(A.16c) 
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c — uqW e 


d^0, 
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— 3a' 


; d0G 
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+ iaRcp (c — uo) 


d0G 

d^ 


+ iaRe, 


dUn 


dz 


(j)G- (A.16d) 


0+ 


Using operator notation to re'write Eq. (A.14)-(A.16) yields 


Ccj) = ici;A40, (A.17) 

which highlights the generalized eigenvalue problem associated with the stability problem. 


3. Numerical method 


We solve this 
location method 


eigenvalue equation numerically by employing a standard Chebyshev col¬ 
li] with an approximation for the streamfunction of the form 


^{z) 



Nl 

^ ^ (^nTn 

n=0 

Na 

y^hnTn 

n=0 



— Sl < Z < 0, 


0 < z < 6 gi 


(A.18) 


where (•) is the Chebyshev polynomial of the first kind. We substitute this trial 
solution into Eq. (A. 17) and evaluate the result at Nl + Ng — 6 Gauss-Lobatto collocation 
points. Together with the eight boundary and interfacial conditions, this yields Nl + Ng + ‘2^ 
linear equations in as many unknowns. In matrix terms, we solve 


Lv = iuMv, (A.19) 

where v = (oq, ..., oat^, 6o) ■ ■ ■) We use a standard linear algebra package (MATLAB®) 

to solve this equation, thereby adjusting the number of collocation points {Nl + 1, Ng + 1) 
until convergence is reached. 
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